{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('/home/luca/GitRepositories/Brancher')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian statistics with Brancher\n",
    "\n",
    "Computers are glorified calculators. Similarly, probabilistic programming is glorified (usually Bayesian) statistics. Therefore, the simplest use of Brancher is as an agile and intuitive toolbox for performing Bayesian inference in few variables.\n",
    "\n",
    "The starting point of all Bayesian models is a probabilistic generative model:\n",
    "\n",
    " $$ p(x|z) p(z)~, $$\n",
    " \n",
    "where $x$ is an observabla variable and $z$ is a latent variable. \n",
    "In this model, $p(x|z)$ (the likelihood) is basically a recipe for generating data $x$ as a stochastic function of an unobservable variable $z$. On the other hand, $p(z)$ (the prior) represents your knowledge of $z$ before observing $x$.  \n",
    "\n",
    "The goal of Baysian inference is to obtain the posterior distribution of $z$ given the model and the observed value of $x$. This posterior is given by the Bayes theorem:\n",
    "\n",
    "$$ p(z|x) = \\frac{p(x| z) p(z)}{p(x)}~. $$\n",
    "\n",
    "It is not trivial to directly use this formula in most real-wold problems (the tricky bit is the innocuous looking term $p(x)$). Brancher deals with this difficulty using an approximate technique called \\amph{stochastic variational inference}. We we deal with the inference later on. For now, let's see how to write generative models in Brancher:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bulding the Beta Binomial Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On of the simplest example of generative model is the Beta-Binomial model. Imagine that you are tring to determine is a coin is fair by tossing it $n = 10$ times. Assume that $\\rho$ is the (unobservable) true probability of the coin landing on head. The probability of observing head $k$ times is given by a Binomial distribution:\n",
    "\n",
    "$$\n",
    "p(k|\\rho) = \\text{Binomial}\\left(n, \\rho \\right)~.\n",
    "$$\n",
    "\n",
    "Now that the likelihood has been specified, we can complete the generative model by specifying a prior. The only strict requirement here is that the probability (density) of $\\rho$ should be cointained in the interval $[0,1]$ since $\\rho$ is a probability. A mathematically convevient choice is the Beta distribution:\n",
    "\n",
    "$$\n",
    "\\rho \\sim \\text{Beta}\\left(a, b \\right)~.\n",
    "$$\n",
    "\n",
    "Let's now write this model in Brancher. The first step is to import the relevant classes of random variables from the standard_variables submodule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from brancher.standard_variables import BinomialVariable as Binomial\n",
    "from brancher.standard_variables import BetaVariable as Beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In brancher, Beta and Binomial are classes that can be used for initializing the random variables that will form the building blocks of our generative model. For example, let's initialize the $\\rho$ variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = Beta(alpha=1., beta=1., name=\"rho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\alpha$ and $\\beta$ are parameters of the Beta distribution. Now that we have $\\rho$, we can initialize $k$ as a binomial variable that depends on $\\rho$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = Binomial(total_count=10, probs=rho, name=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This completes our simple model! We just need to insert the variables we constreucted in a ProbabilisticModel object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from brancher.variables import ProbabilisticModel \n",
    "\n",
    "model = ProbabilisticModel([rho, k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Done, the generative model is now fully encoded as a brancher object!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Observing the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we did not specify that the number of heads $k$ has been oberved. Let's say that we got $8$ heads, we can tell this to brancher by calling the \"observe\" method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k.observe(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting ready for the inference: the variational model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is now to determine the probability of $\\rho$ (yes it is the probability of a probability, a bit confusing...) given the fact that we observed 8 heads. \n",
    "\n",
    "Brancher can solve this problem for you by fitting the parameter of an approximate posterior model. We will refer to this model as the \"variational model\". We need to create a variational model that contains all the non-observed variables in your generative model. Variational models are just probabilistic models and can be created in the very same way. The simplest way to make an appropriate model is to \"copy the prior\". More formally, we create an approximate posterior that follows the same distribution of the prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qrho = Beta(alpha=1., beta=1., name=\"rho\", learnable=True)\n",
    "variational_model = ProbabilisticModel([Qrho])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only difference here is that we set the learnable flag as True. This is important since we need to fit the parameters alpha and beta of Qrho in order to approximate the real posterior. We can now set our variational model as approximate posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.set_posterior_model(variational_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the work is done! Now we can just ask brancher to fit the parameters of the posterior model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/luca/GitRepositories/Brancher/brancher/inference.py:63: UserWarning: The inference method was not specified, using the default reverse KL variational inference\n",
      "  warnings.warn(\"The inference method was not specified, using the default reverse KL variational inference\")\n",
      "100%|██████████| 3000/3000 [00:09<00:00, 309.82it/s]\n"
     ]
    }
   ],
   "source": [
    "from brancher import inference \n",
    "\n",
    "inference.perform_inference(model,\n",
    "                            number_iterations=3000,\n",
    "                            number_samples=10,\n",
    "                            lr=0.1,\n",
    "                            optimizer='Adam')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can wo get the loss and plot it to see if the posterior parameters converged:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f57d123efd0>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8FGX+B/DPN42Q0CF0MLSjKsUoVRQpgqhnL8fdKepxnp797gQ9xXKeeHbPghwK3k9FzoIiUQSkSJESeu8h9ISEEgjpz++Pnd1smd2ZTbbMJJ/365VXdmdnZ59nM/nOM08VpRSIiMg+YqKdACIiCg4DNxGRzTBwExHZDAM3EZHNMHATEdkMAzcRkc0wcBMR2QwDNxGRzTBwExHZTFw4DtqkSROVmpoajkMTEVVLa9euPaGUSjGzb1gCd2pqKjIyMsJxaCKiaklEDpjdl1UlREQ2w8BNRGQzDNxERDbDwE1EZDMM3ERENsPATURkMwzcREQ2Y7nAvWDbcRw/UxjtZBARWZblAve9/83Aje+tiHYyiIgsy3KBGwAOnzof7SQQEVmWJQM3ERH5x8BNRGQzDNxERDZjqcCtlIp2EoiILM9SgZuIiIwxcBMR2Yxh4BaRziKywe3njIg8Eo7EsKaEiMiY4Qo4SqmdAHoBgIjEAjgMYFaY00VERH4EW1UyFMBepZTpJXaIiCi0gg3ctwOYEY6EAABrSoiIjJkO3CKSAOA6AF/4eX2ciGSISEZOTk6o0kdERF6CKXGPArBOKXVc70Wl1BSlVJpSKi0lxdQK80REVAnBBO47EMZqEoADcIiIzDAVuEUkCcBwAF+HNzlERGTEsDsgACilCgA0DnNaiIjIBEuNnGRFCRGRMUsFbiIiMmapwM22SSIiY5YK3EREZIyBm4jIZiwVuMtZV0JEZMhSgbvL03OjnQQiIsuzVOAmIiJjlgrcdRNNjQciIqrRLBW4iYjIGAM3EZHNWCpw5xeWRjsJRESWZ6nATURExhi4iYhshoGbiMhmGLiJiGyGgZuIyGYYuImIbIaBm4jIZhi4iYhshoGbiMhmTAVuEWkgIl+KyA4R2S4i/cOdMCIi0md2Or63AMxVSt0sIgkAksKYJiIiCsAwcItIPQCDAdwFAEqpYgDF4U0WERH5Y6aqpD2AHADTRGS9iEwVkeQwp4uIiPwwE7jjAPQB8L5SqjeAcwDGe+8kIuNEJENEMnJyckKcTCIicjITuA8BOKSUWqU9/xKOQO5BKTVFKZWmlEpLSUkJZRqJiMiNYeBWSh0DcFBEOmubhgLYFtZUERGRX2Z7lTwI4FOtR8k+AGPDlyQiIgrEVOBWSm0AkBbmtBARkQkcOUlEZDMM3ERENsPATURkMwzcREQ2w8BNRGQzDNxERDbDwE1EZDMM3ERENsPATURkMwzcREQ2w8BNRGQzDNxERDbDwE1EZDMM3ERENsPATURkMwzcREQ2Y9nAnXu2KNpJICKyJMsG7psn/xLtJBARWZJlA/f+E+einQQiIkuybOAmIiJ9DNxERDZjapV3EckEkA+gDECpUoorvhMRRYmpwK0ZopQ6EbaUEBGRKawqISKyGbOBWwGYJyJrRWRcOBNERESBma0qGaiUOiIiTQHMF5EdSqmf3XfQAvo4AGjbtm2Ik0lERE6mStxKqSPa72wAswBcqrPPFKVUmlIqLSUlJbSpJCIiF8PALSLJIlLX+RjACABbwp0wIiLSZ6aqpBmAWSLi3P8zpdTcsKaKiIj8MgzcSql9AHpGIC1ERGQCuwMSEdmMrQL3wbwCpI5Px4o9HAdERDWXrQL3msw8AMAXaw9FOSVERNFjq8BNREQWD9xHT5+PdhKIiCzH0oG7/0sLsWhHdrSTQURkKZYO3ACw+fDpaCeBiMhSLB+4iYjIEwM3EZHNMHATEdkMAzcRkc3YMnArpaKdBCKiqLFl4CYiqsksH7hFb5vobSUiqhksH7iJiMgTAzcRkc3YMnCzcZKIajJbBm4ioprMloGbjZNEVJNZPnDrxWhWlRBRTWb5wE1ERJ5MB24RiRWR9SIyJ5wJIiKiwIIpcT8MYHu4EkJEROaYCtwi0hrAaABTw5scc9g4SUQ1mdkS95sA/gagPIxpMY2Nk0RUkxkGbhG5BkC2UmqtwX7jRCRDRDJycnJClkAiIvJkpsQ9EMB1IpIJ4HMAV4rIJ947KaWmKKXSlFJpKSkpIUsgq0WIiDwZBm6l1ASlVGulVCqA2wEsVEr9Nuwp07w+fxcKS8oi9XFERJZn+X7cZeUKU37e57GNpXAiqsnigtlZKbUYwOKwpCSAolLPEjcbJ4moJrN8iZuIiDwxcBMR2YwtArdoC5ixapuIyCaB24lV20RENgvcRERkk8DNKhIiogr2CNzRTgARkYXYInATEVEFBm4iIpuxReB+e+EenC/mfCVERIBNAjcAvLtoT7STQERkCbYJ3O4zBHKSKSKqyWwTuN1jNSeZIqKazEaBm6VsIiLAToE72gkgIrII2wRuADiYdz7aSSAiijr7BG4B3liwK9qpICKKOtsEbgmisiT3bBFSx6dj4Y7jYUwREVF02CZwB2PrkTMAgI+WZUY3IUREYWCbwF2ZTiXsiEJE1ZFtAre3BduOI3V8Og6fikyDZXFpOfafOBeRzyIiCsQwcItIooisFpGNIrJVRJ6LRMK8vb94r+vx6fMlmJlxEACw+dDpKh/79PkSw33+/s1mDHl1MfLOFVf584iIqsJMibsIwJVKqZ4AegEYKSL9wpuswBbtzEFxaTkAoKxc4eHP12NPdr7r9WDGVc7begw9n5uHNZl5AfdbsTcXAHCuqDTo9BIRhZJh4FYOZ7Wn8dpP1MecL9mVAwDYcuQ0vt1wBI9/scn1WjBD4lfucwTsjQdPhTaBRERhYqqOW0RiRWQDgGwA85VSq8KbrODptUOaGSavon8NomruYF4BPlq2P9rJoGrEVOBWSpUppXoBaA3gUhHp4b2PiIwTkQwRycjJyQl1Ov0KVccRzoVibYt2ZNu2X/6Yqavw/JxtOFXA9hEKjaB6lSilTgFYDGCkzmtTlFJpSqm0lJSUECUvcjjjoLWNnb4Gd0/PiHYyKuVMoaPxu5ynGIWImV4lKSLSQHtcG8AwADvCnTCz9ArKwfx/BDMiMxLOF5fhpR+2e8w/Hm4H8wrw7qI9Ibt4DXjpJ0z5ea/xjkRUKWZK3C0ALBKRTQDWwFHHPSe8yTLvbKFnL4+DeQWuyG0mJAeq4/5y7SGkbzpahdQFb+rSffhgyT58tDxydaJjp6/BKz/uxJHThSE53pHThfjn9+au7YUlZeypQxSkOKMdlFKbAPSOQFoq5eNfDrgep286igc+W4e7BqQCCG7kpF4d91++2AgAGH3R6CqlMRjFZY5ujqVlkbuvdq7nGY3qokEvL8SJs8XInBS57zjSrHVPR9WBbUdO6tl0yNGlb9vRM3732X/iHFLHp2OaV4k2K9f8qMizRaVIHZ+ObzccrlxCQyTvXDH+b+UB4x1NikY1/4mzNafBbuOhU1iXdTLayaBqoNoE7g0m+2EPeXUxAOCN+Y4pYg/mFQBwlNzNriR/6KTjPdFewPiRmRvw9DdbsOOY/wsVWcfYaWtw43srop0MqgaqTeD2EEQd94Lt2a7HztGY4bRqXy4yTcx5Yqb0m3euCABQUhqaojJ7RBLZg2Edt92NfnspWjaojf/8Ps30e1LHp+OG3q0M96tMj5TbpqwEAEvW6bJHJJE9VPvAvfXIGdf83EbKlXI10M1ab1x/HapRl1sOn8b5kjJcktrItS2Y0i9HfxLVLNWyqmS1NmFUsKMhe78wH4/9b6PP9vJyhUMnK6aPDXXf72v+vQy3TP4l6PeFOh3VtaqkuLQcJ6M4q2N1H5W7eGc2Jn67JdrJqFGqVeD2Lnca/bus3u87I6BeSfv1+VVf67KguDSoxsxQVFus3p8X1DS0Rp955NR5fLfxSBVTFXkPfLYOvV+YH+1kVFt3TVvj0S2Xwq9aBe5glJUr3PqBuVLunE1VD1avzduFV37cabhfKMtmt37wC24zmUczbpn8Cx6csR5lNhu7PX9bdOc4qd7l7ej6dNWBqPfuioZqFbiDmZr1nMmuf4D/W129EurR0+d1A733CE8jbyzY5XdATGFJGUrKzPWA2Z19NuDr54pKPVYR+nhFJvbl6L/nyOnIrDYUSpGcOiAalu7OwX9/yQy4j1IKM1ZnVcvv4qlZW0wViKqbatU4ucqr6iMagx1u+2AlsvIKMLJ7c8TFBn9ddA/VR08XomWD2j77dHl6LtqnJCM5oep/vodmrPd4PnH2VtStFYfNz11V5WNbwTX/XhbtJITV7z5cDQD4ff9Uv/v8sOUYJny9GZm55zBhVNcIpYzCqVqVuL2dLDBekqwq9AriztJruCsT9uWEZv3LTYcrln5zFvDzDeYOsdNMinsM7jhqgnxtdsJoNtBSaFXrwB1uevFLArxmZNGObI/60J93VX5ec6PbZ6dg0mmjeE0RYnbEMoUWA3clBOrdpVy/g49yY6ev8Xg+/uvN2HksH5OXBD9F6jPfbjXc54Mle3HibJHrudk0Wz1+F5eW4/9+ybRMI2p17Q2YkZmH699dHu1k1EgM3CZU5v/OqHTa/6WfPBoF/bnunWWY9MOOgNUTlS0Jv/SD59SresdZuOO47Zbden/xXjz97VZ8te6Qx3bnd7j7eD7GTF1ZLRvrIsnM+UvhYanAnZQQG9bjrz0Q2sbKYOKld2n26OlCfJlxyM/eFYq0+VMCBeetR86EbaKpu6dn4Pk52zy2Wb3K5KS2RJh3Tx5nup/7bhuW78nV7cdfHTnzbbVFQ6jyLBW4Z90/EE9e3SUsx95/4hxuej/8M7M5S3VKAVm5BcjKLfC/bxChv1wnWjpvwZ+ctRkj31waXEJ101N1VW24XLjjOB6buSEEKfHNj/O583v7/Uersd5Pz6Odx/KrzayL3vkm+7NU4O7cvC7GDe4QlmPPXHMwJMe57F+LXI8D/R8oKAx+ZREGv7LI/z46Me5skf7tu5nq2tMFJRHp8RHoglPVj797ega+NpgnZk/2WdfUunr8BSi97+bfC/UHb1z15s8huRiGUzCjYql6sVTgDqdQBrSiEuPBL94fp/fxeinyt2SZmdJ5z+fn4ZNVWYb7+TN3y7FKvzcU3l9srhF22OtLMOjlRUgdn46tR04bvyEAO3Vt9NbvpZ8Cvq6Uwu7j+RFKjX3syc5H6vh0rNyXG2Cfs8h1a7i3mhoTuLcfq8IJ7FWCu/Ydx6AO57/8yXPFeOmH7Sh1G82YYaY+3UTQcJYenbtuM5jpcNGO7ICv7805i8tfWYTtOqsEvTzX3DqRzrQ8O3srUsene75m6gj63D+/3GSPkK/XHfb7nXgH5XCF6O1HzyB1fDpSx6fjrQW7dfYwrqNYtS8XK/ac0H1NKYVpy/fjVIFnCdto/vgPl+3H8Dd+xoasii5787cdx4SvNxmmpzpbsdcRsAOtJzvs9SW44pXFEUpR8GpM4F5ThYaoQINdlFLo/cJ8fLBkH+ZvO+4KDnd+tNrwuGYCiTP2OOu4F+4IPO+GUQnyx63HcCC3AN9uqPr8K9NXZAb9+YCjkfiswSCfsdPXYOYa47uHD5ftx9VvO6o08s4VY9RbS3EwT7+3w4q9uT4XBL3U/t7E387dqLcqqlTeWOA7IdkJEyW326asxG+mrtJ9bV3WKTz33TaM/2pzUOlar/WxPuhWrfSH/2ZgxurQVBtGUmFJmceFa/Oh07hw4o+VOpbRKersRmo0EC2aakzgLgvTLXGO2z9lUWm57knxsU6AA4KrD3bGG/f3BDv/STA+XeV/trfScoXr3qncUPL8whLc9P4K3P/puoD7LdmVgye+2uxzIfAu4bu78b3l2H70DBZs17+43fnRakxbkekx94x7HM8vLMHhU+erNPDJjF/2+r9F11NU6mj3OFngv0578U7HndbaA3kVk2o5e5NUg0bJ26esRK/n50Npc+ZPWbqv0oHVeU65fy+zNx7BtxsOY/624+jw5PehSHJYGU52ISJtAPwXQHMA5QCmKKXeCnfCQi0cdZmnCkpw6YsV9YwveHWbc5o4W38wTJlS+LXJAKiX/gN5vg10/nKZk1+EeduOmaqfB4AZq7PwjlvDnXvXuewzhdh0SL9u2d/nHz51Hst3n8DQrk0BAJsPmRtxV66AWBOBp6C4FJkBevA47cs56zcwX//ucuwN0VQCgeQEWXdqphvfXdPWIHPSaNz0vmM2yMxJowO2i+w/cQ7tmiQHlY5oco7QvPfjDDRISqjSsVy9bNy2OefsuTWtdZWOHSlmZikqBfC4UmqdiNQFsFZE5iul9KOURYVjEJ33LXBukK38+YUl2OgnAHr7cetxCBwrhTvFCGB2CMl9n6wNqh977tliHD1d6Hr+k1sVTUyAIpz79eXjFZm4rFMT3Dlttav6Ys6Dg0ynAXBUEcWaCFz/mmtuhjjvun33C2IkgjbgCBLX9WwZkc9y94Nb4/OQVxdjWNemGNPvAgzp3LRSx6tqWSi/sATZ+UXokFLH9Ht+0tpwqvL9PfedI3SFeoGL9xbvQU5+ESZe2z2kx9VjWFWilDqqlFqnPc4HsB2A8YKMFmOV4c/uPllpvgfIX77YiMe/2OixuHFJmW+e/P0zebeQG52z8V4zG36wZJ/rsV7dtlP65oq684mzt2LB9uMedc7O2frMTgBm9u9mtnpD5ysLm8KSMkxduk/3tX/M2Wa6ATZUTp/3/M4XbM/G2Glr/OxdISu3IKjeO4UlZaa6Kt72wUoMfW2J6eMGUlJWjmdnbzXVnhCM9E1HTRd4/jV3J6Ytzwzp5/sTVB23iKQC6A1AvxWFoq60vBxHdIYiBzvFbGyM/8geKHA/OtN36beqMFuq23fCt7Ss+17vniba01UBuoadq2Rd6ruL9uAf6dt1X5u6bD+26fTsCSTYMB+qEZODX1mE0W+bb9N4aMZ6XO423sGfYPMfyIJtxzF9RaarNG1GQbHv39X7nHngs3WmBu7p/c+Fk+kJnUWkDoCvADyilPL5xkVkHIBxANC2bduQJZCCs3xPLgZMWuizPS5AINYTqjaBQIcx6toI6I8YrQrvqikFhYN5Bbhtykrd/dM3HcUDn63DnAcHoUer+kF9Vn6IGo9dd0dBfhUqjI2TSilk6bSxAMC8KKw4VKrdvXy38Qjq1IrD8j0ncLaoFF/e1x/t/VTF+Et/Zfxk0A031EwVw0QkHo6g/alS6mu9fZRSU5RSaUqptJSUlFCmkUIgUL20nmDu4h+buQEdnvw+6GDv7MYXSLh6Azkt35MbsES1ZJfjH7KqA32qwvmXW52Zh9Tx6VBK6c6z4ttvvXLf3Yo9J5A6Ph0HvQKb+ziFr9cd9lmLVSmFBSaDtnvPmlBUY/57YUX/+Rmrs5CVV4C8c8UB7w7PnDd/YS0vVzgd5vn9g2EYuMVRg/8hgO1KqdfDnySKBKNRisGUdL9efxhl5Ur3H7CqpU5VXvmqCrOBKzvff71oOK8bO4/lGw6i0fP5moOm1kutKHEbX7RX7DmB9E1H8f3mo5iZ4ejn7V232/GpH1yP1x/0rfctKi3Hvf/N8NhWUFyK699d7jPvyx3/qbjDMbsMXyC7jusvmBGowLLt6BmfVbL8FRRu+eAX9Hx+nqm7xEgwU+IeCOB3AK4UkQ3az9VhThdF2aGTwdfZleoE7owDVZuBr0wpXB5gvpdA9BpvK6u4tByTl+x1lTqzcgv8rs3pZBQvH/9iIybO3uKxTe+uxbuhz9/n+kyzEPjjXUrKyvGbqavwwGfrPPrX6134nEFWL77p/f1X7c/DhoOn8NL3/kflOt+3+dBp3PT+ChSWlKGwpMxVEHjAoM9/IIFK3ABw43ted1t+vjTnRWx3tjWmEDDTq2SZUkqUUhcppXppP9bvoU4e9us03oVaqOujncc8cbZykymFsu/+5CX7MOmHHXjgs3UoKSvH4FcW4coAPSJ2Hsv3qWrQM2P1QY9RpIt3+vaO+ZNX4PJXQHX//svKlVvjZGDed0qBVnF6+hvHhUbvmy0NUHJ27v/Mt1tw8QvzPV5zpvuZ2Vuw9sBJbD1yBl2enouHPnf0rU7f7H9oeqgZnTHOu5ev1h5CRmb0pgWuMSMna7LJS/bifAQWDdArcZnlL8hW5WLw6jzf4efBcn66M7j+uPU4XnM7rr+qjqve/Nmj62Yg7lVBZwqN61H9fSfuAe7v32xx1c8blfyD+YqdjXB67yn2Ctzni8tcF4EirRT9318O+Ix38O4W6bw4BJpLJBjvLNSbP0af0fnmbON//IuNuHmycXVVuDBw1wCTfjA3eVRV6fVLNtsVrd0E/Zu48qpXf1aJsxHNvQ+0+1Jy3+uUBoOts3Xfv6xc4dnZW3Eg1/8dkr/gssVt4ecZq7NcVUVGf4FgLo4xAszdctRndSHAt2qq6zNzcZfWT3zV/jx0eXqu7jF7Pe9ZAjfTTTCYtvZX5+0yPc2AUdnjTd1JxICcM4W628PFVoH7sz/0jXYSKIDndfrQVrWK5o35VS81V4XR8lx69cCd3BrxzHCvqthy+Aymr8jEnz9b73d/f4HWX/w1GkDifTz3xkzvnhSxIrjvk3W6dxqBLjaV9d1G/cnQgp0kzWyfcaOL2J7ssx4NlEoppI5Px9t+5nUPF9P9uK0g2C5tFFl6CyAcq2JJxNnDoTqb7RaEyrRbjM2HT6OguBQn8n3r9/2VCv2FnDMGPXv8HW/tgZN44ivPKWAD9VB5rYpVU3pHDvUc8SPeMBipaeLmw70b6/Ez0Zmz25Il7p6t9Qc6MGyT1ZwvrnpdzmtudxXu9eLbj+bjrum+U8z6Gypf2eYAf2tvfroqy6f6I1A1UFXWdPXXxrEkxDM16nUbHPRyxYC1YNtUorXKvSUD9zcPDNTdHupJYYiq6ulvtxjvFAT3qpnCkjLdueALikPb0PwHr77XswIsHReoz3tV/PXLTViX5TtjpNG87aHg3vU12ItfVe8oK8uSgdtfgGbcJqsJ5+Rl7rMzupvtp943HN0xI+XLtb6NnaFktmtoZUebRpolA7c/jNtUkyzdHVw1gdFgEzK2bLf+8nFWY6/AzSI31SChWF6OHN5ZtAeHThoPiDoX4mqocLFZ4K543LJ+YvQSQkS2cqqgxHC5PDuxVeBOiI3BVd2bRTsZRGRD/pbbsyPLB+67B7ZzPY4RwdPXdDN8T+342HAmiYgoqiwbuJvWrYW/j+6KZ67thi7N6wJwVJW4T1V53+UddN+7/pnhaNWgdqSSSkQUUZYN3KufGoZ7L2sPAKillaC9R06OH9VF972J8bFokBQf3gQSEUWJZQO3u/fG9MHDQzvhV83MrwY9YVTXMKaIiMhXckJkqmltMVdJqwa18ejwX5na99Fhjv0GdWoSziQREfmoFaH2NVsEbnetG9bGXQNSMaav/oLEAzs2jnCKiIgiy3aBW0Tw7HXd/b7eu23DsH5+83qJUZufgIgIsEkddzBiYwKPrnR/OaVuraCPP+ehQUG/h0KvXqLtyhxUA0RqbHe1C9xG4mIrsvzjI4ODfn+jpIRKf/aIbsEPHorGgKOmlbigRVrdxOB6DdWtxUBP4RepWTlsH7hHdm9uuE+ftg1cj90H5zRKTkCHlGSPfS/r1ASv3dLT5xjP/7o7Zt0/ADE6Jfq9//Rd9P7G3q18tt0/pKNhWr29cH0Pn22NkxPw0JXBH8usIZ2bhu3YoVIrzvyp2zg5ASufHBrG1BA5RSZyG579IvKRiGSLSGgnHg6Ryb+7GJmTRgfc56s/DXA9fv7XnvXjM/7QD++P6eN6/n/39EUnrdthx6YV3Q+v69nSb/25XvVMS50BQLFul+PhJkvfifGx+HxcP49tr97SE4+N6Gzq/e76tmvkely3Vhw2Thyhu9/jV5nrwRMKbRpVbqBUfKz5wP3DI5chuVYcrr7Q+CJv1lu39wrZsYiCZebsnw5gZJjTEVYigsxJo5E5aTQaelV1NK2XiFEXtvDYFhfj+Fpa1E9E/drBD+R5f0wfPDysk046HL+7tahn6rq8+dkRqJcYjzpet/lDulSuRHxjH7e7AIHfvDWtG94JvNyrf1LqVK5aJj7OfMnG+Rnef/tAhnROCfj6r3v53lGFwgWNk0zt17ZREr64r39Y0hBKN1/c2u9rdpjs09m92GoMA7dS6mcA+msbVVNdW9TFxGu74c3bApeqmtSphdVPOW7BZ90/AGMHpmLLc1dh1IUtPEqEcx4chL+N7Owa+VmuPKdrb9tI/5/VWY9rdIK7VwU5De3SFKO1C1JCbAxevKEHhnZ1K+VrCWhWryJwfnpvX2x6Vr8U7o/3HcxdA1Jdj/u39+2aOW5we7x+a8X3qjdtvccFxksLbVZIf+uPznnQs/G4dcParumAgwncL990kcfz3/Rti8/u7YuebXy/ayPu3wkAJMTFoKGfkb0vXn+h3+PMur/izjFG4JoKIlTc70yD1ShZ/7vt1NT/oLmLWgf/XZoRH2v+ivDaLT3Rr30jv6/rFcCswPZ13MHq4OdEevWWnljwmKOxUkQwdmA7NK5TK+DKGRl/H+YqnfZu2xATr+3uUzoGgB6t6uP+Kzq6qlTcVyp56MqOmPfoYNzYpxVuTdMvnRgtkvz1/Y6l3q7onIKb+jiOcX3vVnjjtl745oGB2PXiKIzpe4Fu2v55gyNQ9G/fGAM7NkE9r0a/j+5KC/jZ7qXznf8YieRaFW0I4y53TFlwqVsVTdtGSUjW0tG6YW2P9oTH/AyyeuO2nnhGm1zsklTHsUQE6To9fHq08lyvtENKxd/7zwbtAkPd7mTidKpiBnRsgm/uH4B9Om0agCP/zu/f3fVe7R0xAix8/ArdYwSqOnKvfjuQVxDytVoCffaOF0Zi2RND/L7ur5dPoFPX/c7LqLrTiPv4jcdHdMbT13RD2gUN0T4lGQkBqtVuurg1Prmnr+v58G7N8Bs/Y0SsJGRN7SIyDsA4AGjb1hoZ//juS3Hea2L0Vg0cA3gGdfQcWenvlm7a2Evw8YoDPgEtmKu60wWNk9CkTgImjOqKT1dlAQC6t6qPxPhYvH5rL5SVK9wzqL3fQekAAAANOUlEQVTPye7v5I+LEZRqS2c5T/yT54rRKDkeI3s0R3xsDHq5lRAT42Ox8ZkR6Pn8PFc9vnsQ1ZMYVxGIlz0xBENfW4Ki0ooFY93TFiOCpATH8R4Z1glJWkOw3sVvwWOXo0mdBDRwKwU7S/+xIvjj5e3xwZJ9AIAberdGWblC3cQ4XNS6AWZvPIIOKcno3rIiSG9+doTrwvTW7b2Qe7YY7VOS0eeCinaJxACj2nq3bYAP77oEqePTXWnQIyJ+/x614mLxys0X4at1gZfhihFBw+QETP7txZi75Si+0VkwYUCHxlixN9f1/MJW9dGsXiJaNaiNw6fOY/rYS/0uGuy044WR6PL03ID71IqLQVFpud9Cw+w/D8TX6w6jVlwMWjf0X42j4Cg4LN6Zg5duvBAbsk5hZsZB3c/r07Yh1madxNgB7bA28yTGDXZc4P/3x/4YM3WlzwLFZrifYndc0hb1k+JxzyDHzKIvpm/Df5bu9/veuNgYbH9+JF6euwN/vaozvtt4BJ8F8dlJCbGY/edB+HbDYZ+LdLiELHArpaYAmAIAaWlpEV+47aLW9X3m2738V/r1lIEG8Hi7+IJGuPgCz1upDc8M1+1dYiQxPhYZfx8OAK7A7X6U2BhBZ53bX/cS98RrK6a1XfrEEGSf8Vy8tWFyAp4a7X/q2/pJ8XhiZBdX42jfdo3wr5suwjU9W+ju7/xD9m/fGK0bJrmC1sCOjbF8T65H+gXAPYPaoaxc4b7LO2DjIcfir+7/VFdodccdde58KmZ+dMw1s3JfHjYedBwjNkZwS1obAMD0sZegbzvPahj37oFm65+Hd2uGPw5uj5sn/+LT5pBgotfK5N9ejA+X7cOazIrVzfXOC+8L14NXOm6/R/ZojuKyco/AnRAXg4WPX44mdWuhsLgMl/7zJwBwNVAv+ssVKCotQ93EeBSXBl5hPjE+FgseG4w/fbIOu7N9VzcHgL+M6IytR07jr1d1ca23mBAXg7kPX4a4mBi0bZwUsErjqu7N8OPW41AKeP3WXjiYV4CebRrgsLYAb604z4vltLGXYECHikLTh3dd4np8abtG2P3i1VBK4ZsNh/HozI1o06g2lvxlCNo/+T0uSW2IQR1TkL75CPbmnPNY7/OG3q2wYm8uXvh1d9T3qoZ6dPiv8OGy/fC+zr14Q0WPrdoJsa64MLBjcNNlzHt0MFo3TMLjlegwUFnVpnPr//7YH+cisCI0AI9SYiDf/XlQSGYpdNbN3jUgFWPd5idvUb82WtQPvlfGn66omA5XRHDrJW387tujVX0kJcTioaGOYHNR6wZYvT/PVRKPc7vziBFBYnyMa9/uLeuhXZNkjB/VBc3rJ6KguEy31PaP63ugfUqy65/dWWr/z+8u1l1V/Aq37or/vOFCj+oZI+/+pg9+3pWDmRkHESviujB5L4tXOyEWSQmxAVdUH9mjOUb2aO4qpTu9fmtPPPa/ja7n7vHCu0rgup4tcUGjJBQUlyHvXLHH37NeYjxm/MERsJ13RglxMa6LSkJcDDInjUZBcSnWZ53CmKmrAAD92jdytbF0bFrX406jc7O62Hk8H9f2bIkxfduib7tGrrxnayOC6yXGo32KfpVi1xb1sP3oGdfzJ6/u6gjcUGiUnOCq635gSEfECHDHpW0xcfZW1/7Ohv9ARAQ39G6NvHMlGNGtGWJiBD88fBlaN6yNuonxeGhoR/y49Rju+8Sxok3DpHjcktbGdWH3lpQQh6VPXImBkxa6tt1xaRuM6XuB7v5t/LQ5ufvqTwNw0/srACDgnUi4GAZuEZkB4AoATUTkEICJSqkPw52wYCXGxwa8FY6GC1vX9/tav/aNsGD7cbQ10YsgpW4tLP3bEFfDXCSs/fswFJaWo37teGx7vqJT0dQ707D7+Flc0DgJ//l5Hy7/VUUQ9a5CSEqIw6K/XGH4Wb/t5/gHKiwpw45j+a551pvWS0TTeoHzHGx95OiLWmBYt6YoLVf461WdkZXnWIfQmfTl4690FQDWPT3cVdVwZ/9Uv8f87s+D0LhOxcVc7x//t/3a4pOVWbrvD9Tg2b+D8dw7SQlxHqXEz8d59jZx7/P+6PBOuO+TdWicnIB+Oo3HRj4f1w9ZuQW49p1lAADRvjnv2rDaCbG6XVbrBjHi1VnVATguGE4ighStbeXuge3w16uMS7ru8/NveGa4bnuPu0/v7Yv1WY47qfSHBiH7TBE6pNTBR8v3Y/qKTHRMqYPbL2mDDdodYaQZfotKqTsikZCa5p5B7TCyR3PTV2szpYBQauynm169xHhcrNUbT7jaMXXuDb1bYdb6w1VezDkxPtbUCkdVVSsuFq/d6mgUdVYPXK31wHH/B4/Tqj3GDW6vW4Xl5H2B7tWmAW7s3Qrzth3H2aJSKAX84/oL8Y8APUZCYVjXZvh5l+/K8G/f0RsDJi1Ep6Z1cEXnphg3uD3udQuKwahfOx7dW7oHUcfvAG34GNO3ratqsLnBhdisPm0b4D+/T8PAjo1R2+RUqsO7NUPX5nVN3TEP7NjEdTHs3rI+urd0bH/2uu6uKpVJXj2PIkkC9ZqorLS0NJWRkRHy45I1lZaVo6CkzKcB1y4KiktROz5W98JTXq4g4luVYsbBvAJMXboPz1zb3XAOnVA4X1yG/MISw7sUf3Lyi3DJiwvQtlESfv6b/x4kSim0m/A9gIoG0JduvBB3XOr/7kcphcKSctNBtiYSkbVKqcDduJz7MnATEeAIru8t3ovRF7ZAapNk4zdQSAUTuKtN4yQRVY2I4IFKzKdDkVfjBuAQEdkdAzcRkc0wcBMR2QwDNxGRzTBwExHZDAM3EZHNMHATEdkMAzcRkc2EZeSkiOQAOFDJtzcBcCKEyYmm6pKX6pIPgHmxouqSD6BqeblAKRV4zTxNWAJ3VYhIhtlhn1ZXXfJSXfIBMC9WVF3yAUQuL6wqISKyGQZuIiKbsWLgnhLtBIRQdclLdckHwLxYUXXJBxChvFiujpuIiAKzYombiIgCsEzgFpGRIrJTRPaIyPhop0ePiHwkItkissVtWyMRmS8iu7XfDbXtIiJva/nZJCJ93N5zp7b/bhG5Mwr5aCMii0Rku4hsFZGHbZyXRBFZLSIbtbw8p21vJyKrtHTNFJEEbXst7fke7fVUt2NN0LbvFJGrIp0Xt3TEish6EZmjPbdlXkQkU0Q2i8gGEcnQttnxHGsgIl+KyA7tf6Z/1POhlIr6D4BYAHsBtAeQAGAjgG7RTpdOOgcD6ANgi9u2fwEYrz0eD+Bl7fHVAH6AYx3afgBWadsbAdin/W6oPW4Y4Xy0ANBHe1wXwC4A3WyaFwFQR3scD2CVlsb/Abhd2z4ZwJ+0x/cDmKw9vh3ATO1xN+28qwWgnXY+xkbpPHsMwGcA5mjPbZkXAJkAmnhts+M59jGAe7XHCQAaRDsfET8p/Xwx/QH86PZ8AoAJ0U6Xn7SmwjNw7wTQQnvcAsBO7fEHAO7w3g/AHQA+cNvusV+U8vQtgOF2zwuAJADrAPSFYxBEnPf5BeBHAP21x3HafuJ9zrnvF+E8tAbwE4ArAczR0mbXvGTCN3Db6hwDUA/AfmjtgVbJh1WqSloBOOj2/JC2zQ6aKaWOAoD2u6m23V+eLJVX7fa6NxwlVVvmRata2AAgG8B8OEqYp5RSpTrpcqVZe/00gMawSF4AvAngbwDKteeNYd+8KADzRGStiIzTttntHGsPIAfANK36aqqIJCPK+bBK4NZbAtvu3V385ckyeRWROgC+AvCIUupMoF11tlkmL0qpMqVULzhKq5cC6Kq3m/bbsnkRkWsAZCul1rpv1tnV8nnRDFRK9QEwCsADIjI4wL5WzUscHNWj7yulegM4B0fViD8RyYdVAvchAG3cnrcGcCRKaQnWcRFpAQDa72xtu788WSKvIhIPR9D+VCn1tbbZlnlxUkqdArAYjrrFBiLiXAzbPV2uNGuv1weQB2vkZSCA60QkE8DncFSXvAl75gVKqSPa72wAs+C4qNrtHDsE4JBSapX2/Es4AnlU82GVwL0GQCet9TwBjoaW2VFOk1mzAThbiO+Eo77Yuf33WitzPwCntVuqHwGMEJGGWkv0CG1bxIiIAPgQwHal1OtuL9kxLyki0kB7XBvAMADbASwCcLO2m3denHm8GcBC5ah0nA3gdq2nRjsAnQCsjkwuHJRSE5RSrZVSqXD8DyxUSo2BDfMiIskiUtf5GI5zYwtsdo4ppY4BOCginbVNQwFsi3o+It1gEaAR4Go4ejfsBfBUtNPjJ40zABwFUALHFfQeOOoUfwKwW/vdSNtXALyr5WczgDS349wNYI/2MzYK+RgEx23aJgAbtJ+rbZqXiwCs1/KyBcAz2vb2cASrPQC+AFBL256oPd+jvd7e7VhPaXncCWBUlM+1K1DRq8R2edHSvFH72er8n7bpOdYLQIZ2jn0DR6+QqOaDIyeJiGzGKlUlRERkEgM3EZHNMHATEdkMAzcRkc0wcBMR2QwDNxGRzTBwExHZDAM3EZHN/D/tbCazESDf3QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "loss_list = model.diagnostics[\"loss curve\"]\n",
    "plt.plot(loss_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Everything seems to be fine. We can finally plot the resulting (approximate) posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAACsCAYAAACXSStGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGIZJREFUeJztnXmUXHWVxz+3el+zbyQkYckyAYPBqCgR2Q0Mg4rMQEZUGIRhRkUFHVEZUByHMzqjZwAZ4cgioKyKRJYBhGRYA4kJSUggMZAQsqez9b7WnT/eq6rXlerqqk6/fvVe3c85ffJ7v/d7r253+tv3t937E1XFMIzwEAvaAMMw8sNEaxghw0RrGCHDRGsYIcNEaxghw0RrGCHDN9GKyJ0isktE3uzjvojITSKyQURWicjxftliGFHCT097NzA/y/2zgGnu1+XA//hoi2FEBt9Eq6ovAHuzNPk0cI86LAGGi8gEv+wxjKgQ5Jh2IvC+53qLW3cQInK5iCwTkWXHHHOMAvYVri9jEAlStJKhLuN/sKrerqpzVXVuVVWVz2YZRmETpGi3AId7ricB2wKyxTBCQ5CiXQh80Z1FPgE4oKrbA7THMEJBqV8vFpH7gZOB0SKyBbgeKANQ1V8CTwJnAxuAVuASv2wpKro7YN97MHoaSKYRiBF2fBOtqi7o574CX/Hr84uSrja44wzYsRpmngMX3GfCjSC2IypKvP2EI1iAtx+HLUuDtcfwBRNtlNi2ovf1uieDscPwFRNtlNi+svf15iXB2GH4iok2Suzd2Pt6+0qI9wRji+EbJtqooArNO3vXdbVCw/pg7DF8w0QbFdr2Qbzr4Pr0ca4Reky0UaFpR+Z6E23kMNFGhWYTbbFgoo0KLQ2p8vjZqfKO1dDTPfT2GL5hoo0KbftT5fqJUD3aKXe3w54Nwdhk+IKJNiq0H0iVK2ph5BGp6x2rht4ewzdMtFGh3eNpy2tg5JGpaxNtpPAtYMAYYnqJthbqDktdJ/YjG5HAPG1UaEsTbS9PmzEhphFSTLRRwTumLa+F2rEQK3OuWxt63zdCjYk2KqSPaSUGdZ7klnveGXqbDF/wVbQiMl9E1rkJya/JcH+yiCwSkRVuwvKz/bQn0qR7WoB6z7h277tDa4/hG36eMFAC/AInKfksYIGIzEprdi3wkKrOAS4EbvXLnsjT5hVtjfNv3bhU3YH3MaKBn572I8AGVX1XVTuBB3ASlHtRoN4tD8OyMQ6MeBw6GlPXCdFWjUrV9bU32Qgdfoo2l2TkPwAuchO/PQl8LdOLvMnKd+/e7Yet4abjAMmU0WVVECtxytVe0Vqiy6jgp2hzSUa+ALhbVSfhZGa8V0QOssmbrHzMmDE+mBpyMo1nAapHpsrmaSODn6LNJRn5pcBDAKr6KlAJjPbRpmjSljZznKCXaM3TRgU/RbsUmCYiR4hIOc5E08K0NpuB0wBE5K9wRGv933zpy9NWpXlatWN1ooCfp+Z1A18FngbewpklXiMiN4jIuW6zq4HLRGQlcD9wsZsP2ciH9C2MCcqqoKzaKfd0OtktjNDj695jVX0SZ4LJW3edp7wWONFPG4qCvjwtON62q9UpN+3o3WU2QontiIoCfY1pASrrU+XWPUNjj+ErJtookB5L66XCI9q2bGd8G2HBRBsF0vcdezFPGzlMtFEg25jW62lbzdNGARNtFMg2pq2oS5VNtJHARBsFcva01j2OAibaKNCeIcInQaVNREUNE20U6GtzBZinjSAm2ihgE1FFhYk27HS1OwnJAWKlUFrR+75NREUOE23YSR/PSlpEZCJfFDhxtz0ZTtYzQoWJNuxkG8+CExDvrbeggdBjog072WaOE9hkVKQw0Yad9CTlmajwiNnb3gglJtqwk4un9YrZkpaHHhNt2OlvTJte326eNuwEmqzcbfN3IrJWRNaIyG/9tCeSZIvwSdZ7J6JMtGHHt8wVnmTlZ+AkeVsqIgvdbBWJNtOA7wInquo+ERnrlz2RxSvC9FjaBOZpI0XQycovA36hqvsAVHWXj/ZEk2y7oRLYRFSkCDpZ+XRguoi8LCJLRGR+phdZsvIs5D2mtYmosBN0svJSYBpwMk7i8l+JyPCDHrJk5X2T9+yxedqwE3Sy8i3AY6rapaobgXU4IjZyJZd1WpuIihRBJyv/A3AKgIiMxuku25mM+ZDTjijztFEi6GTlTwN7RGQtsAj4tqraPrt8yHdMa5429ASdrFyBq9wvI1/icWjPcMRlOjamjRS2IyrMdDSS8YjLdMqqUuF5Xa3Q3Tkk5gWBiKiI3Ou5LhWR3SLyeJ7v2eQO2QqOnEQrIpNE5FH3m98pIr8TkUl+G2f0Qy5rtODE2Hq9cLSXfVqAY0Wkyr0+A9gaoD2DTq6e9i6cSaQJOGutf3TrjCDJZQtj8n5RdZGfAv7aLS/AOdwNABEZKSJ/EJFV7t6A2W79KBF5RkRWiMhteJYsReQiEXldRN4Qkdvc3X6Bkatox6jqXara7X7dDdiCadDk6mnT70d/MuoB4EIRqQRmA6957v0QWKGqs4HvAfe49dcDL6nqHBwHNRmSR7BegLPV9oNAD/D5Ifku+iDXiagGEbmI1F+sBYDN8gZNLmu0me5H3NOq6ioRmYrze/pk2u15wOfcds+7HnYYcBJwnlv/hIgkUnycBnwIZ+88QBUQ6HbbXEX7D8AtwM9xZj5eceuMIMlljTZBRdFtZVwI/CfObrtRnvpsO/UynY0swK9V9buDat0hkFP3WFU3q+q5qjpGVceq6mdU9T2/jTP6IZc12uR9b9BAUeSJuhO4QVVXp9W/gNu9FZGTgQZVbUyrPwsY4bZ/Djg/EYHmjomn+G9+3+TkaUVkDE5EzlTvM6pq3jZI8vG03vsdjX23iwiqugX47wy3fgDcJSKrgFbgS279D4H7RWQ58H/AZvc9a0XkWuAZEYkBXcBXgMCcVq7d48eAF4E/4QzEjUIgl1jaBGXFseSjqgf9IFR1MbDYLe/l4BBR3J14Z3qqvum59yDw4CCbOmByFW21qn7HV0uM/MnL01Z7nou+p40yuS75PC4iZ/tqiZE/+Yxpy4qrexxlsnpaEWnCmT2rci6lA6dPLzhbh+uzPW/4jHnaoiSraFW1TpzFqT+r6vFDZJORK/ms05qnjQz9do/dSJxXReTDQ2CPkQ95LfmYp40KuU5EnQJcISKbcDZkJ7rHs/0yzOgH1fxmj4tsySfK5DoRdRZwJHAq8DfAOe6/Wckl77Hb7nw3pGpujvYYXW3Q0+GUY2VQUpG9fa8ln+ISrYj0uJv93xSRh0Wkuo92T2bKUVZo5ORpB7L7KZe8x267OuBKem/qNvqjPc3Lph9xmU5ZVarc2QTxnr7jb6NHm7vZHxH5DXAF8LPETXfeRlQ15xUSzzPxwTa2P4LOewzwI+AnQLuPtkQP71bE8rq+2yWIlfQWbkfT4NsUDl4EjhaRqSLylojcCiwHDvcGvovIVa5nflNEvuHWHfRMEN9AoHmPRWQOcLiqZs0qYHmPM5DPeDZBkc8gi0gpzlAvsR95BnCPqs7x9iZF5EPAJcBHgROAy9zf1T6fGUoCy3vs7uP8OXB1fy+yvMcZ6OVpcxRt8c4gV4nIG8AynD3Fd7j176nqkgzt5wGPqmqLqjYDvwc+0c8zQ4afid36y3tcBxwLLHbjFMcDC0XkXFVd5qNd0SB9TJsLRbL/OAPJMW0C93eupY/22SYI+npmyAgs77GqHlDV0ao6VVWnAksAE2yuHKqnLcLucR68AHxGRKpFpAb4LM5YuCAIOu+xMVB6jWlzmIiCol72yQdVXQ7cDbyOs6rxK1VdEahRHgLNe5xWf7KftkQO87Q500e43iac4Zm3bqqn/DM8y0J9PRMElvc4rAxkTFs8aVQjjYk2rAzE0xb5kk9UMNGGlYGMaYt3ySdSmGjDStveVNk8bVFhog0rLQ2pcmWOe9zLbfY4Cphow0hnK3Q2O+VYaf9ZKxIU6exx1DDRhpFWr5cd1n+ET4IiXafNNTSvn3d8Y4DP3SAip+f7XDZMtGGk2RM0kWvXGIo5EL5NVT+oqscCnTihefnyDSAv0YpIiapep6p/yueZ/tqYaMNIi0e0VfmI1jt7XLTrtC8CR0Of4Xc1IvKEiKx06y8QkSuBw4BFIrLIbXemiLwqIstd713r1m8SketE5CXgb0XkbhE53713mnsq32oRuVNEKjI90983YKINIy2e858qh+X+XKk3EL7ZCYQvIryheVnC7+YD21T1ONcz/6+q3oQT7HKKqp7ixtxeC5zuJjxcBlzl+ah2VZ2nqg94PrsSZ2vkBar6AZzdiP+U7Zm+MNGGkZYBdo8PCoQvmi5yptC8vsLvVgOni8h/iMgnVDVTl+QEYBbwsvveLwHe830ynUYwA9ioquvd61/jnNSX7ZmM+Lr32PCJgSz3JCircfJLgTMZVTUie/tokCk0L+Psnaqud73w2cCNIvKMqt6Q1kyAZ1V1QR+flyl8r7/ZwpxD/szThpHGralyvqKzZZ8EGcPvROQwoFVV78M5KjOR77sJJwYcnDDSE0UkMTauFpHp/Xze28DUxDPAF3AO+sob87RhZL8ni09tnpk8inTZJx1VXS4id+OE34EbficinwJ+KiJxnNM0EuPO24GnRGS7O669GOeUvUQazGuB9fSBqraLyCXAw+7Yeinwy4HYbqINIwc8oq0Zm9+zRehpM4XmufWZwu+exokBT297M3Cz5/p54KAE/t7wPvf6Yk/5OWBO2iMHPdMf1j0OG11tqYkoiUH1qOzt0/HuUy5iTxtmfBVtf8nK3XWytSKySkSek4BP2A4FBzzj2erR+ecuLis+Txs1fBOtJ1n5WTjT4wtEZFZasxXAXPd4kUdw8h8b2TiwOVWuGUBmSguEDz2BJitX1UWq2upeLsHJ2GhkY+/GVLk2z/EsFOWYNmoEmqw8jUuBpzLdsGTlHnavS5WHDyDBvc0eh57AkpX3aihyETAX+Gmm+5as3EODR7TDJuf/vHna0BNksnIA3LCl7wOfVNUOH+0ZMAfautjd1E5bZ5zh1WWMra+gojT/w6samjtYu62Rtdsb2XGgnab2bjq6e6ivKmN4VRmTR1YzY3wd08fVUVPRx3+N19MOO1RPa2PaMOKnaJPJyoGtOMnK/97bwN2gfRswX1V3HfyK4OjqifPA65t5ZPlWVm3Zj3r6CDGBSSOqOXJMDUeOrnX+dct1laW0dHbT0NTJX3Y1sWFXM2u2NbJm2wF2Nub+N2nyyGpmjq9j5oR6Zo6vY8b4OiZWdlLZtN01ohTqxuf/jVmeqNDjm2hVtVtEEsnKS4A7E8nKgWWquhCnO1yLs0sEYLOqBp7I/M2tB7jqoTdYv7M54/24wua9rWze28ridf6MsRPvf2btzmTdibHV/KbcKW+SSdz4bAc1ZUJ1GZSXCKrO+KMnriTOX4wBFaXC2Gph5sgSPlxRTTJkwLrHoSTQZOWqOqgR/YPByxsauPyeZbR0psLWYgJj6yqpKI3R2N7F/tauzIPzfigviTF5VDVTR1Uzvr6K6ooSykpitHZ209jWxbb97by/r5Xt+9vp0YM/YY5sSJZf6jiKpzd1523D1Jiy2BW+edpwYtsYPfz5vb1cctdSOnscP1VRGuO8ORM5deY4aitTP6rO7jg7GtvZvr+NbQfa2X6gje3uv909SkVZCTXlJUwYVslhw6uYMqqGqaOqmTCsipJY/6lhunribNvflvS2m/e0snV/G3M630m2WRGfNqDvcXc8tSOqp6WBWDyOxGxjXJgw0brsONDOFfctTwp2ZE0518yfyeEjD84wUl4aY/LIaiZnuDcYlJXEmDKqhimjPJNGqsx9eBO4w+J5HziK6aXQ1gPt3dAVd6brS8RJGZUIPFOF9h7Y1Qpv7oHNzVU0ayW10k5JvJMbH13CNed9LHGKnBECTLRAPK5cef8Kdjc5iqitKOW6c2Yxrr4yYMtSVDS/T1mHk+u4p7Sa6YdPYPoAdPbmHti7fAS1OBNazy9bTUnNSP5l/szBNNfwEesXAQ8sfZ/XNzmCiAl8/bRpBSVYgLqG1KFtbcOOcoIFBsCxo2DE8FTg/DjZx62L3+HBpZuzPGUUEkUv2l2N7dz41FvJ63OPO4xjJ+aRd2mIqG1YmSy31h+dpWX/9FSmAufH4pwJdP3CNazb0XRI7zWGhqIX7Q8fX0tTuzMLO76+ks/OKcztz7W7vZ720ETbXZES7fRqZ1mrvSvOV367nPau4kr2FkaKWrTPv72TJ1ZtT15fOu8IyksL70ciPR3U7Ev1BtqGHXVI7+suT3WPT5sUp8L9njfsauanT6/r6zGjQCi839AhoqWjm3/9w5rk9UnTRhdktxigZu9bxOKdAHRUjaOnvP6Q3tfl8bQjexr4wgmpMOY7XtrIK+80ZHrMKBCKVrQ/f3Y9W/c7WQlrK0r5/AmFG39ft3tZsnyoXWOArupUSF9V4yZOnTmWDx6e8r7ffngVje1dh/w5hj8UpWhXbznAnS+n4lK/+LEp1FeWBWhRdup3Lk2WW0cc+tJMR/VhyXJl0yZE41x+0pHUukEKW/e3ccMf1x7y5xj+UHSi7eju4VsPryTu7hI8duIw5h09OlijsqFx6nalRNsy/NBFGy+rpssd18binVS0bGFEdTmXzjsi2eaRP2/h6TU7DvmzjMGn6ER703N/Yd1OZ2mjojTGl+cdUdC7gar2/4WyTufU9+6yOjprDuvnidzwvqeq8V0ATjhyFCd6/oB97/eraWguyGjJoqaoRLt43S5uXZzav7vgI5MLbhNFOsO3v5gst46Ymfuxlv3QUTMhWa7Zm+oKX/zxqYyscSIK9rR0cs3vVqMZgheM4Cga0W5saOHK+1ck42JnTajnjFnjgjWqP1QZvfGx5GXT6INS5g6Y1mGpgIO6XamJrtqKUv7xpCOT1396aye3PL8Bo3AoCtFu2NXEBbe9SqO7iWJkTTlfO/VoYgXcLQYY//bd1O51lqXisTKaxs4dtHd7J7Tqdy2DeCrMb/ak4Zzp+YP2X8+u5/YX3jGPWyBEWrSqymNvbOW8W19hlxsMUF4S45unT2d4dXk/TwfLqI1/ZOqyf0te75t4Kj1lGRPlD4iuyjF0VYwEoKS7xRGuh4tOmNJr3frfn3ybqx5ayb6WzkGzwRgYQScrrxCRB937r4nI1EP9zHhc2djQwn1L3uOcm1/i6w+8kfSwFaUxvnPWTI4eO3i//H4wfOsijn7lW4gbat867Gh2TrtgcD9EhKYxxycvx6+7t9ftspIYV58xnRnj6pJ1j67Yykk/WcSPHl/La+/uoaUj/yB849ARv7o8brLy9cAZOEnelgILVHWtp80/A7NV9QoRuRD4rKpm/e2cO3euLlvmeIWHlr7PYyu30tLRQ0tHN62dPexu7qCzO37Qc2PrKvjaqdMKVrDDty5mzDuPUNGyvVdET3vNRDZ9+PpB9bIJKhs3ctRr309eN4+aTWfVGHYf9Tn2Tp4POAH/d728kcXrD06rIwIjq8sZUVNOXWUpZbEYZaXCzQuOT05mJZoOuvFFjJ/xtMlk5QAikkhW7l21/zTwA7f8CHCLiIjm+Jfk/X2tvLxhT9Y25SUxPnXMOC46YQrV5YUbPlzfvp3R7/XKzENX1Ri2z/sxJVWjyT/3Y//Ea46jcdInqd/inLhYu2cVAB0TP05nIhtkBVx95gw+ftQo7lnyHlv2tSWfV3VmmPekdZnjNvb1FT897fk4WRa/7F5/Afioqn7V0+ZNt80W9/odt01D2rsuBy53L2cAg72rfTRgG24d/PhZNKjq/EF+Z9Hip+vJJVl5TgnNVfV2nPNBfUFElqnq4E3Nhhj7WRQ+fk5E5ZKsPNnGPWh3GLDXR5sMI/T4KdpksnIRKcdJVr4wrc1C4Etu+Xzg+VzHs4ZRrASdrPwO4F4R2YDjYS/0y55+8K3rHULsZ1Hg+DYRZRiGP0R6R5RhRBETrWGEjKIWbX/bLIsJEblTRHa5a+dGAVO0onW3Wf4COAuYBSwQkVnBWhUodwO2ASIEFK1o8WyzVNVOILHNsihR1RewNfJQUMyinQi877ne4tYZRkFTzKLNaQulYRQaxSzaXLZZGkbBUcyizWWbpWEUHEUrWlXtBhLbLN8CHlLVNdmfii4icj/wKjBDRLaIyKVB22RkxrYxGkbIKFpPaxhhxURrGCHDRGsYIcNEaxghw0RrGCHDRBsgItIctA1G+DDRBoQ452vaz9/IG/ulGUJEZKqIvCUitwLLgSoR+bGIrBSRJSIyzm03RUSeE5FV7r+Tg7XcKCRMtEPPDOAeVU2cW7lEVY8DXgAuc+tucdvMBn4D3DT0ZhqFiu2IGkLcA8YWqeoR7nUHUKmqKiIXAGeo6pdFpAGYoKpdIlIGbFfV0X2+2CgqzNMOPS2ecpcnz3MPfae0tb+sRhITbWHyCqkc0J8HXgrQFqPAKNxj5IqbK4E7ReTbwG7gkoDtMQoIG9MaRsiw7rFhhAwTrWGEDBOtYYQME61hhAwTrWGEDBOtYYQME61hhIz/BxfHYxsjv8AuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 260.375x180 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from brancher.visualizations import plot_posterior\n",
    "\n",
    "plot_posterior(model, variables=[\"rho\"]) #TODO: to be fixed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Pytorch",
   "language": "python",
   "name": "pytorch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
